function kummer2

%  plot airy functions

% clear all previous variables and plots
clear *
clf

n=200;
x=linspace(-5,5,n);

%alpha=1; beta=1;
%xt=1; yt=-0.75;
%axis([-5 5 -1 1])
%loc='Northwest';

%alpha=-0.5; beta=2;
%xt=-2; yt=-8;
%axis([-5 5 -10 10])
%loc='North';

alpha=1; beta=-0.5;
xt=-1.2; yt=-2.5;
axis([-5 5 -3 3])
loc='East';

%alpha=0.5; beta=2;
%xt=-1.7; yt=-0.8;
%axis([-5 5 -1 1.01])
%loc='Northwest';

% get(gcf)
set(gcf,'Position', [824 372 515 280]);

for ix=1:n
	z=-0.5*alpha*x(ix)^2;
	a1=0.5*beta/alpha;  b1=0.5;
	a2=0.5*(alpha+beta)/alpha;  b2=1.5;
	y1(ix)=M(a1,b1,z);
	y2(ix)=x(ix)*M(a2,b2,z);
end;

hold on
box on
grid on
plot(x,y1,'-','Linewidth',1)
plot(x,y2,'--r','Linewidth',1)

say=['\alpha = ',num2str(alpha),', \beta = ',num2str(beta)]; 
text(xt,yt,say,'FontSize',18,'FontWeight','bold')

xlabel('x-axis','FontSize',16,'FontWeight','bold')
%ylabel('Kummer Solutions','FontSize',16,'FontWeight','bold')

%set(gca,'xtick',[0 1 2]);
%set(gca,'ytick',[0 0.5]);
%set(gca,'XTickLabel',{'0';'t_M';'2t_M'})
%set(gca,'YTickLabel',{'0';'x_M'})

set(gca,'FontSize',14);
legend(' y_1',' y_2','Location',loc);
set(findobj(gcf,'tag','legend'),'FontSize',16); 


function g=M(a0,b0,x0)

% this code will not work if b0 is a negative integer
%if (b0==round(b0)) & (b0<0)
%	error('The code doesn't work if b is a negative integer');
%	return;
%end

% use identity:  M(a,b,x)=exp(x)*M(b-a,b,-x) 
if x0< 0
	x=-x0; a=b0-a0; b=b0;
else
	x=x0; a=a0; b=b0;
end

xmax=20;

% evaluate M(a,b,x) for x ge 0

%  use exact expressions for M if possible
if (x==0) | (a==0) 
	g=1;
elseif a==b
	g=exp(x);
elseif (a==1) & (b==2)
	g=(exp(x)-1)/x;

%  using the complete poly isn't the smartest idea if a is a very large 
%  negative integer but its done anyway 
elseif (a==round(a)) & (a<0)
	N=-a;
	g=1;
	cs=1;
	for in=1:N
		cs=cs*(a+in-1)/((b+in-1)*in);
		g=g+cs*x^in;
	end

%  otherwise, for small x use series definition
elseif x<xmax
	N=800;
	g=1;
	cs=1;
	for in=1:N
		cs=cs*(a+in-1)*x/((b+in-1)*in);
		g=g+cs;
		if abs(cs/g)<1.0e-12 
			break
		end 
	end

%  and use asym approx for large x
else
	N=40;
	g0=1;
	cs=1;
	for in=1:N
		cs=cs*(b-a+in-1)*(1-a+in-1)/(in*x);
		g0=g0+cs;
		if abs(cs/g0)<1.0e-12 
			break
		end 
	end
	q1=gamma(b)/gamma(a);
	g=x^(a-b)*exp(x)*q1*g0;
end

%  adjust answer if x0<0

if (x0<0)
	g=g*exp(x0);
end





